from vwellnopr import *
from vwellpr import *
from datetime import *

class IntEnsBombeoCC(VWellNoPR): 
  """
  Clase de base para el calculos de parametros de pozo en flujo no permanente

  Funciona para todos los acuiferos
  """
  def __init__(self, tipoAc=None,  Q=None, r=None, t=None,  s=None,  alpha=None,  rdB=None, h0=None):
    """
    Inicializador de clase
  
    La decripcion larga
  
    Parameters
    -----------------
    tipoAc : string
      L = Libre, C = Confinado, SC = semiconfinado.
    Q : float
      Caudal de extraccion durante el ensayo en m3/d.      
    r : double
      Distancia desde centro de bombeo hasta el punto de observacion en metros.
    t : array float
      Valores de tiempo de las observaciones en minutos      
    s : array float
      Valores de descensos observados en metros .
    alpha : array float
      alpha = (S coeficiente de almacenamiento)/(T transmisividad).
      Por defecto se toman valores de S entre 0.5 y 1x10-5 y de T de 1 a 2000 m2/d
      obteniendo valores de alpha entre 5x10-1 a 1x10-9
    rdB : array float
      rdb = (r distancia desde el bombeo a la observacion)/(B).
      Por defecto se toman valores de S entre 0.001 y 2.5
    h0 : float
      h0 = espesor saturado inicial del acuifero
    """
    PozoNoPerm.__init__(self)
    if tipoAc:      
      self.tipoAc = tipoAc    
    if Q:
      self.q = Q
    if r:
      self.r = r
    if t:
      self.t = t
    if s:
      self.s = s
    if alpha:
      self.alpha = alpha
    else:
#      if self.tipoAc == 'L':
#        self.alpha = self.generaRangoLog(0.000000001, 5.0, 30)
#      else:
#        self.alpha = self.generaRangoLog(0.000000001, 0.5, 30)
      self.alpha = self.generaRangoLog(0.000000001, 5.0, 35)
    if rdB:
      self.rdB = rdB      
    else:
      self.rdB = self.generaRangoRDB()
    if h0:
      self.h0 = h0  

  def setTipoAc(self, tipoacuifero):
    self.tipoAc = tipoacuifero
    
  def setCaudal(self, caudal):
    self.q = caudal
    
  def setDistancia(self, distancia):
    self.r = distancia
  
  def setData(self, t, s):
    self.t = t
    self.s = s


  def generaRangoRDB(self):
      rdB = ['0.001','0.002','0.003','0.004','0.005','0.006','0.007','0.008','0.009','0.010','0.020','0.030','0.040','0.050','0.060','0.070','0.080','0.090','0.100','0.200','0.300','0.400','0.500','0.600','0.700','0.800','0.900','1.000','1.100','1.200','1.300','1.400','1.500','1.600','1.700','1.800','1.900','2.000','2.100','2.200','2.300','2.400','2.500']
      return ( rdB )
  
  def generaRangoLog(self, li, lf, n):
    """
    Genera intervalos equilog 
  
    La decripcion larga
  
    Parameters
    -----------------
    li : float
      Limite inferior
    lf : float
      Limite superior      
    n : float
      Nuero de intervalos en el rango
    """    
    # from math import *
    b = []
    kk = ( lf/li )**(1.0/n)
    for i in range(0, n+1):
      b.append( li*kk**i )
    return ( b )    

#  def generaRangoRDB(self):
#    b = []
#    for i in range(0, 4):
#      for j in range(1, 10):
#        b.append( j/10.0**(3-i) )
#    return ( b )

  def generaRango(self, nf=4, mf=10):
    """
    Genera intervalos equilog 
  
    La decripcion larga
  
    Parameters
    -----------------
    li : float
      Limite inferior
    lf : float
      Limite superior      
    n : float
      Nuero de intervalos en el rango
    """    
    # from math import *
    b = []
    for i in range(0, nf):
      for j in range(1, mf):
        b.append( j/10.0**(3-i) )
    return ( b )
    
  def obtenerVTransAcC(self, a):
    """
    Obtiene el vector de transmisividades para ese valor de a = S/T    

    Parameters
    -----------------
    a : float
      Valor de a = S/T
    """
    trans = []    
    u = []
    cont = 0
    for t in self.t: 
      u.append( self.r*self.r*a/(4*t) )
      # TODO: usar funcion
      trans.append( self.q * self.Wu( u[cont] )[0] / ( 4.0 * self.pi()*self.s[cont] )  )
      cont = cont + 1
    return ( trans )

  def obtenerMTransAcC(self):    
    """
    Obtiene la matriz de transmisividades para todos los valores de a = S/T    
    """    
    trans = []
    for a in self.alpha:      
      trans.append( self.obtenerVTransAcC( a ) )      
    return ( trans )

  def obtenerTransAcC(self):    
    return ( self.obtenerMTransAcC()     )

  def obtenerVTransAcL(self, a):
    # self.h0*self.q * self.Wu( u[cont] )[0] / ( 2.0 * self.pi()*( self.h0**2 - ( self.h0 - self.s[cont]  )**2 )  )
    # H0*k = T
    """
    Obtiene el vector de transmisividades para ese valor de a = S/T    

    Parameters
    -----------------
    a : float
      Valor de a = S/T
    """
    trans = []    
    u = []
    cont = 0
    for t in self.t:
      # u = r2*S/(4*h0*k*t) 
      u.append( self.r*self.r*a/(4*t) )
      # TODO: usar funcion
      trans.append( self.h0*self.q * self.Wu( u[cont] )[0] / ( 2.0 * self.pi()*( self.h0**2 - ( self.h0 - self.s[cont]  )**2 )  ) )
      cont = cont + 1
    return ( trans )

  def obtenerMTransAcL(self):    
    """
    Obtiene la matriz de transmisividades para todos los valores de a = S/T    
    """    
    trans = []
    for a in self.alpha:      
      trans.append( self.obtenerVTransAcL( a ) )      
    return ( trans )

  def obtenerTransAcL(self):    
    return ( self.obtenerMTransAcL()     )

#  def  obtenerCondAcL(self):
#    """
#    a = S/k
#    """
#    k = []
#    for a in self.alpha:
#      u = []
#      kalpha = []
#      cont = 0
#      for t in self.t: 
#        u.append( self.r*self.r*a/(4*self.h0*t) )
#        # TODO: usar funcion
#        kalpha.append( self.q * self.Wu( u[cont] )[0] / ( 2.0 * self.pi()*( self.h0**2 - ( self.h0 - self.s[cont]  )**2 )  ) )
#        cont = cont + 1
#      k.append( kalpha )      
#    return ( k )

  def obtenerVTransAcSC(self, a, rdb):
    """
    Obtiene el vector T para un a = S/T y rdb = r/B
    """
    trans = []
    cont = 0
    for t in self.t:
      u = self.r*self.r*a/(4*t)
      trans.append( self.q * self.Wua( u , rdb  )[0] / ( 4.0 * self.pi()*self.s[cont] ) )
      cont = cont + 1
    return ( trans )


  def obtenerTransAcSC(self):  
    trans = []
    for rdb in self.rdB:
      transrdb = []
      for a in self.alpha:
        u = []
        talpha = []
        cont = 0
        for t in self.t: 
          u.append( self.r*self.r*a/(4*t) )
          # TODO: usar funcion
          talpha.append( self.q * self.Wua( u[cont] , rdb  )[0] / ( 4.0 * self.pi()*self.s[cont] ) )
          cont = cont + 1
        transrdb.append( talpha )
      trans.append( transrdb )
    return ( trans )
  
  def distRelativaV(self, v, cmM = 1.8, minV = 0.0001):
    """
    Funcion que encuentra la distancia relativa entre los elementos de un vector    
    
    Parameters
    -----------------
    v : array of float
      El vector
    cmM: float
      Cociente entre el maximo y minimo del vector
    minV: float
      Valor minimo de un elemento del vector    
    """
    su = 0
    l = len(v)
    vmax = max(v)
    vmin = min (v)
    # los puedo juntar 
    if  vmin < minV:
      su = 2000
    else:
      # Acoto el cociente entre max y min para hacer mas eficiente el algoritmo
      if  vmax > cmM*vmin:  
        su = 2000
      else:
        for i in range(0, l):
          for j in range(i+1, l):
            su = su + ( v[i]-v[j] )**2 / ( v[i]+v[j] )**2
    return ( su )  

  def distRelativaM(self, m, cmM = 1.8, minT = 0.001):
    """
    Funcion que encuentra la distancia relativa entre los elementos de un vector    
    
    Parameters
    -----------------
    m : array of array of float
      El vector
    cmM: float
      Cociente entre el maximo y minimo del vector
    minV: float
      Valor minimo de un elemento del vector    
    """    
    d = []
    for v in m:
      su = self.distRelativaV(v, cmM, minT)
      d.append( su )
    return ( d )
   
  def distRelativa(self, m, cmM = 1.8, minT = 0.001):
    """
    Ver distanciaRelativaM    
    """    
    return ( self.distRelativaM(m, cmM, minT) )

  def media(self, v):
    """
    Funcion que encuentra la media de los elementos de un vector    
    
    Parameters
    -----------------
    v : array of float
      El vector    
    """    
    return ( sum( v ) / len( v ) )

  def calcParametrosAcC(self, alpha=None):
    if alpha != None:
      self.alpha = alpha
    ai = self.alpha[ 0 ]
    af = self.alpha[ (len( self.alpha ) - 1) ]
    Ti = self.media ( self.obtenerVTransAcC( ai ) )
    Tf = self.media ( self.obtenerVTransAcC( af ) )
    rango = self.alpha
    #print "(ai,af,Ti,Tf)=(%E,%E,%E,%E)" % (ai, af, Ti, Tf)
    iter = 0
    while max( [Tf, Ti] ) >= 1.01*min( [Tf, Ti] ): # 1% de error
      distV = []
      cont = 0
      transM = []
      for a in rango:
        iter = iter + 1 
        transV = self.obtenerVTransAcC(a)
        transM.append( transV )
        distV.append ( self.distRelativaV( transV ) )
        # TODO: Mal
        if cont == 1 and distV[cont-1] < distV[cont]:
            ai = 0.1*rango[cont-1]
            af = rango[cont]
            Ti = self.media ( transM [ cont - 1 ] )
            Tf = self.media ( transM [ cont ] )
            rango = self.generaRangoLog(ai, af, 30) 
           #print "Ac. Libre: Break - Distancias (cont=%d,a= %f): (%f,%f)" % ( cont, a,  distV[cont-1], distV[cont] )
            break # siguiente paso del loop while        
        if cont > 1:
          if distV[cont-2] > distV[cont-1] and distV[cont-1] < distV[cont]:
            ai = rango[cont-2]
            af = rango[cont]
            Ti = self.media ( transM [ cont - 2 ] )
            Tf = self.media ( transM [ cont ] )
            rango = self.generaRangoLog(ai, af, 30) 
            #print "Distancias (cont=%d): (%f,%f,%f)" % ( cont, distV[cont-2], distV[cont-1], distV[cont] )
            break # siguiente paso del loop while
          else:
            #print "Distancias (cont=%d): (%f,%f,%f)" % ( cont, distV[cont-2], distV[cont-1], distV[cont] )
            pass  
        cont = cont + 1
   #print "Se hicieron %d iteraciones" % iter
   # La matriz de transmisividades
    self.mT =  self.obtenerTransAcC()      
    # El vector de distancias relativas
    self.vD =  self.distRelativa( self.mT )    
    self.selAlpha =  ai
    #TODO: Ver pesos
    # Hallo la media del vector T seleccionado
    self.trans = Ti
    # Con T y alpha calculo el S
    self.almac = self.selAlpha * self.trans 
    #print "Sali de calcParametrosAcC"
    self.sCalc = self.descAcC(self.t, self.q, self.r, self.trans, self.almac)
    # calculo los descensos con estos valores
    self.ttCalc = self.generaRangoLog(self.t[0], self.t[len(self.t) -1], 30) 
    self.ssCalc = self.descAcC(self.ttCalc, self.q,self.r, self.trans, self.almac)

  def calcParametrosAcL(self, alpha=None):
    di = datetime.now()
    if alpha != None:
      self.alpha = alpha
    ai = self.alpha[ 0 ]
    af = self.alpha[ (len( self.alpha ) - 1) ]
    Ti = self.media ( self.obtenerVTransAcL( ai ) )
    Tf = self.media ( self.obtenerVTransAcL( af ) )
    rango = self.alpha
    iter = 0
    while max( [Tf, Ti] ) >= 1.01*min( [Tf, Ti] ): # 1% de error
     #print "Ac. Libre: **** (ai,af,Ti,Tf)=(%E,%E,%E,%E)" % (ai, af, Ti, Tf)
      distV = []
      cont = 0
      transM = []
      for a in rango:
        iter = iter + 1 
        transV = self.obtenerVTransAcL(a)
        transM.append( transV )
        distV.append ( self.distRelativaV( transV ) )
        # TODO: Mal
        if cont == 1 and distV[cont-1] < distV[cont]:
            ai = 0.5*rango[cont-1]
            af = rango[cont]
            Ti = self.media ( transM [ cont - 1 ] )
            Tf = self.media ( transM [ cont ] )
            rango = self.generaRangoLog(ai, af, 30) 
           #print "Ac. Libre: Break - Distancias (cont=%d,a= %f): (%f,%f)" % ( cont, a,  distV[cont-1], distV[cont] )
            break # siguiente paso del loop while          
        if cont > 1:
          if distV[cont-2] > distV[cont-1] and distV[cont-1] < distV[cont]:
            ai = rango[cont-2]
            af = rango[cont]
            Ti = self.media ( transM [ cont - 2 ] )
            Tf = self.media ( transM [ cont ] )
            rango = self.generaRangoLog(ai, af, 30) 
           #print "Ac. Libre: Break - Distancias (cont=%d,a= %f): (%f,%f,%f)" % ( cont, a,  distV[cont-2], distV[cont-1], distV[cont] )
            break # siguiente paso del loop while
          else:
           #print "Ac. Libre: Sigo - Distancias (cont=%d,a= %f): (%f,%f,%f)" % ( cont, a, distV[cont-2], distV[cont-1], distV[cont] )
            pass  
        cont = cont + 1
   #print "Se hicieron %d iteraciones" % iter
   # La matriz de transmisividades
    self.mT =  self.obtenerTransAcL()      
    # El vector de distancias relativas
    self.vD =  self.distRelativa( self.mT )    
    self.selAlpha =  ai
    #TODO: Ver pesos
    # Hallo la media del vector T seleccionado
    self.trans = Ti
    self.cond = self.trans/self.h0
    # Con T y alpha calculo el S
    self.almac = self.selAlpha * self.trans 
    df = datetime.now()
    d = df -di
   #print "calcParametrosAcL: Demore %d segundos y %d microsegundos" % ( d.seconds, d.microseconds )    
    # Que volores tengo de descenso para esos parametros en los tiempos observados
    self.sCalc = self.descAcL(self.t, self.q, self.r, self.cond, self.almac, self.h0)
    # calculo los descensos para un intervalo de tiempo 
    # aceptable. OJO: s = Ho-H
    self.ttCalc = self.generaRangoLog(self.t[0], self.t[len(self.t) -1], 30) 
    self.ssCalc = self.descAcL(self.ttCalc, self.q, self.r, self.cond, self.almac, self.h0)

  def calcParametrosAcSC(self, alpha=None, rdB=None):
#    if alpha != None:
#      self.alpha = alpha
#    if rdB != None:
#      self.rdB = rdB
    iter = 0
    # Para almacenar [T]
    matRDBTrans = []
    # Para almacenar T de cada [T]
    vectorRDBT  = []  
    # Para almacenar alpha de cada T
    vectRDBAlpha = []
    # Para distancia de cada [T]
    vectRDBDist = []
    for rdB in self.rdB:
      distV = []
      cont = 0
      transM = []
      for a in self.alpha:
        iter = iter + 1 
        transV = self.obtenerVTransAcSC(a, rdB)
        transM.append( transV )
        distV.append ( self.distRelativaV( transV ) )
        cont = cont + 1      
      dmina = min( distV )      
      jj = distV.index( dmina )
      # El valor de la distancia minima
      vectRDBDist.append( dmina )      
      # El alpha correspondiente
      vectRDBAlpha.append( self.alpha[jj] )      
      # El vector con la distancia minima
      matRDBTrans.append( transM[jj] )
    self.vD = vectRDBDist
    # El minimo de las distancias
    self.distMin = min( self.vD )
    # Tomo numeral para ver cual es el alpha adecuado
    ii = self.vD.index( self.distMin )     
    # Encuentro el alpha 
    self.selrdB =  self.rdB[ii]    
    # Encuentro el alpha 
    self.selAlpha =  vectRDBAlpha[ii]
    # TODO: Ver pesos
    # Hallo la media del vector T seleccionado
    self.trans = sum(matRDBTrans[ii])/len(matRDBTrans[ii])
    # Con T y alpha calculo el S
    self.almac = self.selAlpha * self.trans 
    # Calculo los descensos con estos valores
    self.sCalc = self.descAcSC(self.r, self.t, self.q, self.trans, self.almac, float( self.selrdB ), 'calc')
    # calculo los descensos para un intervalo de tiempo 
    # aceptable. OJO: s = Ho-H
    self.ttCalc = self.generaRangoLog(self.t[0], self.t[len(self.t) -1], 30) 
    self.ssCalc = self.descAcSC(self.r, self.ttCalc, self.q, self.trans, self.almac, float( self.selrdB), 'calc')

#  def calcParametrosAcSC(self, alpha=None, rdB=None):
#    if alpha != None:
#      self.alpha = alpha
#    if rdB != None:
#      self.rdB = rdB      
#    # TODO: Re-hacer esto
#    # La matriz de transmisividades
#    self.mT =  self.obtenerTransAcSC()
#   #print "mT = %s\n" % self.mT
#    # La matriz de distancias relativas
#    mDist = []
#    # El vector de distancias minimas
#    vDistMin =  []
#    # El vector con indice de distancias minimas
#    vIndDistMin =  []
#    for i in range(0, len( self.rdB ) ):
#      vD = self.distRelativa( self.mT[i] )
#      minvD = min (vD)
#      ii = vD.index( minvD ) 
#      # Un vector distancia para cada rdB
#      mDist.append( vD )
#      # Un vector con las distancias minimas
#      vDistMin.append( minvD )
#      # Un vector con el indice alpha correspondiendo a esa distancia minima
#      vIndDistMin.append ( ii ) 
#    self.mD = mDist
#   #print "mD = %s\n" % self.mD
#    self.vD = vDistMin
#   #print "vD = %s\n" % self.vD
#    self.vIndexAlpha = vIndDistMin
#    # El minimo de las distancias
#    self.distMin = min( self.vD )
#    # Tomo numeral para ver cual es el alpha adecuado
#    ii = self.vD.index( self.distMin ) 
#    # Encuentro el alpha 
#    self.selrdB =  self.rdB[ii]    
#    # Encuentro el alpha 
#    jj = vIndDistMin[ii]
#    self.selAlpha =  self.alpha[jj]
#    # TODO: Ver pesos
#    # Hallo la media del vector T seleccionado
#    self.trans = sum( self.mT[ii] [jj] ) / len( self.mT[ii][jj] ) 
#    # Con T y alpha calculo el S
#    self.almac = self.selAlpha * self.trans 
#    # calculo los descensos con estos valores
#    self.sCalc = self.descAcSC(self.r, self.t, self.q, self.trans, self.almac, self.selrdB)
#    # calculo los descensos para un intervalo de tiempo 
#    # aceptable. OJO: s = Ho-H
#    self.ttCalc = self.generaRangoLog(self.t[0], self.t[len(self.t) -1], 30) 
#    self.ssCalc = self.descAcSC(self.r, self.ttCalc, self.q, self.trans, self.almac, self.selrdB)

#  def calcParametrosAcL(self, alpha=None):
#    if alpha != None:
#      self.alpha = alpha
#    # La matriz de transmisividades
#    self.mk =  self.obtenerCondAcL()      
#    # El vector de distancias relativas
#    self.vD =  self.distRelativa( self.mk, minT=0.001)
#    # El minimo de las distancias
#    self.distMin = min( self.vD )
#    #Tomo numeral para ver cual es el alpha adecuado
#    ii = self.vD.index( self.distMin ) 
#    #Encuentro el alpha 
#    self.selAlpha =  self.alpha[ii]
#    #TODO: Ver pesos
#    # Hallo la media del vector T seleccionado
#    self.cond = sum( self.mk[ii] ) / len( self.mk[ii] ) 
#    self.trans = self.cond*self.h0
#    # Con T y alpha calculo el S
#    self.almac = self.selAlpha * self.cond
#    # Que volores tengo de descenso para esos parametros en los tiempos observados
#    self.sCalc = self.descAcL(self.t, self.q, self.r, self.cond, self.almac, self.h0)
#    # calculo los descensos para un intervalo de tiempo 
#    # aceptable. OJO: s = Ho-H
#    self.ttCalc = self.generaRangoLog(self.t[0], self.t[len(self.t) -1], 30) 
#    self.ssCalc = self.descAcL(self.ttCalc, self.q, self.r, self.cond, self.almac, self.h0)

  def calcParametros(self):
    if self.tipoAc == 'L':
      # Calculo
      self.calcParametrosAcL()
    if self.tipoAc == 'C':
      # Calculo
      self.calcParametrosAcC()
    if self.tipoAc == 'SC':
      # Calculo
      self.calcParametrosAcSC()
